Spontaneously ordered motion of self-propelled particles 
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We study a biologically inspired, inherently non- 
equilibrium model consisting of self-propelled particles. In 
the model, particles move on a plane with a velocity of con- 
stant magnitude; they locally interact with their neighbors 
by choosing at each time step a velocity direction equal to 
the average direction of their neighbors. Thus, in the limit of 
vanishing velocities the model becomes analogous to a Monte- 
Carlo realization of the classical XY ferromagnet. We show 
by large-scale numerical simulations that, unlike in the equi- 
librium XY model, a long-range ordered phase characterized 
by non- vanishing net flow <j> emerges in this system in a phase 
space domain bordered by a critical line along which the fluc- 
tuations of the order parameter diverge. The corresponding 
phase diagram as a function of two parameters, the amplitude 
of noise n and the average density of the particles g is calcu- 
lated and is found to have the form n c (g) ~ g 1 ^ 2 . We also 
find that <f) scales as a function of the external bias h (field or 
"wind") according to a power law <f> ~ h 9 . In the ordered 
phase the system shows long-range correlated fluctuations and 
1// noise. 



I. INTRODUCTION 

Recently there has been an increasing interest in the 
studies of far-from-equilibrium systems typical in our 
natural and social environment. Concepts originated 
from the physics of phase transitions in equilibrium sys- 
tems [1] such as collective behavior, scale invariance 
and renormalization have been shown to be useful in 
the understanding of various non-equilibrium systems as 
well. Simple algorithmic models have been helpful in 
the extraction of the basic properties of various far-from- 
equilibrium phenomena, like diffusion limited growth [2] , 
self-organized criticality [3] or surface roughening [4]. 
Motion and related transport phenomena represent a fur- 
ther characteristic aspect of non-equilibrium processes. 
Indeed, the transport in various driven systems, such as 
traffic models [5], molecular motors [6] and other self- 
propelled systems [7-12] have been the subject of recent 
studies. 

Self-propulsion is an essential feature of most living 
systems. Moreover, the motion of the organisms is usu- 
ally controlled not only by some external fields, but also 
by interactions with other organisms in their neighbor- 
hood. In Ref. [7], a simple model was introduced cap- 
turing these features with a view toward modeling the 



collective motion of large groups of organisms [13-16] 
such as schools of fish, herds of quadrupeds, flocks of 
birds, or groups of migrating bacteria [17-21]. The aim 
of this paper is to further investigate the various interest- 
ing phenomena exhibited by this novel non-equilibrium 
model. 

The model consists of particles moving on a plane and 
characterized by their (off-lattice) location afj and veloc- 
ity Vi pointing in the direction $j — 0(#j), where the 
function gives the angle between its argument vector 
and a selected direction (e.g., horizontal coordinate axis). 
The magnitude of the velocity is fixed to vq to account 
for the self-propelled nature of the particles. A simple lo- 
cal interaction is defined in the model: at each time step 
a given particle assumes the average direction of motion 
of the particles in its local neighborhood S(i) with some 
uncertainty, as described by 



i (i + At) = <0(t)>s W +6 



(1) 



where the noise £ is a random variable with a uniform 
distribution in the interval [—77/2,77/2] and the local av- 
erage direction of motion ($)s(i) is defined as 



J2 g i 



(2) 



v.. 



The local surrounding of the ith particle S(i) will be 
specified in Sec. II. The locations of the particles are 
updated in each time step as 



Xi(t + At) = Xi(t) + Vi(t)At. 



(3) 



This model is a transport related, non-equilibrium ana- 
log of the ferromagnetic models, with the important dif- 
ference that it is inherently dynamic: the elementary 
event is the motion of a particle at each time step and 
a change in the direction of motion. The analogy is as 
follows: the Hamiltonian tending to align the spins in the 
same direction in the case of equilibrium ferromagnets is 
replaced by the rule of aligning the direction of motion 
of particles. The amplitude of the random perturbations 
is in analogy with the temperature [7]. Indeed, if vq = 0, 
the model is similar to the Monte-Carlo simulations of 
diluted XY ferromagnets [22]. 

The reported long-range order in the above [7] and 
the closely related models [10,11] is surprising, because 
in the case of equilibrium systems possessing continuous 
rotational symmetry the ordered phase is destroyed at 
finite temperatures [23]. A recent dynamic renormaliza- 
tion group treatment of the problem [9] by Tu and Toner 
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has also led to the conclusion of the existence of an or- 
dered phase in two dimensions. Thus, the question of how 
the ordered phase emerges due to the non-equilibrium na- 
ture of the model is of considerable theoretical interest as 
well. In Section II we study the kinetic phase transition 
leading to the symmetry-broken state, which is charac- 
terized in Section III. In Section IV. we investigate the 
effect of an external field applied to the system. 



II. KINETIC PHASE TRANSITION 

We studied the behavior of the model defined through 
Eqs. (l)-(3) by performing large-scale Monte-Carlo sim- 
ulations as a function of two control parameters: the 
density of particles g and amplitude of the noise n. We 
applied random initial conditions and periodic bound- 
ary conditions. The calculations were performed on a 
Connection Machine 5 parallel computer, with typically 
TV = 10 4 — 10 5 particles. 

The interaction range of S(i) was defined in two differ- 
ent ways: (i) as a circle of radius R, or (ii) by considering 
a square lattice on the plane built up from lattice cells 
of length R, and assuming that a given particle inter- 
acts with all the particles located in the same lattice cell 
and in the eight neighboring cells (sec Fig. 1). The exis- 
tence of the long-range order, and the critical exponents, 
turned out to be robust against changing these details 
[(i) or (ii)] of the interaction. Here we present results 
obtained by definition (ii), as this latter choice of S(i) 
increases the speed of the simulations by a considerable 
amount. 

A natural dimensionless parameter is C = v n At/R, 
and the behavior of the model on t>o, At and R depends 
only through their combination given in this expression 
for C. Thus, we can work with At = 1 and R = 1. 
The results presented in the following were obtained with 
v a = 0.1, the role of this velocity is discussed in Sec. V. 

For the statistical characterization of the model, a well- 
suited order parameter is the magnitude of the average 
momentum of the system: 



1 

TV 



(4) 



This measure of the net flow is non-zero in the ordered 
phase, and vanishes (for an infinite system) in the disor- 
dered phase. 

We start the simulations from a disordered system 
(random positions and orientations), thus (f)(t = 0) « 0. 
After some relaxation time a steady state emerges indi- 
cated, e.g., by the convergence of the cumulative average 
(1/t) Jq (f>(t)dt. Here we focus on the statistical proper- 
ties of the steady state only, and do not deal with the 
entire relaxation process. In the vicinity of the critical 
regime to reach the stationary behavior takes more than 
10 4 Monte-Carlo steps for a typical simulation with 10 5 



particles in a 100 x 100 system. In such cases we run 
the simulations for « 10 5 time steps, which takes about 
4 hours CPU time on the CM5. 

The stationary values of <j) that we obtained as a time 
average are plotted in Fig. 2a vs n for g = 2 and various 
system sizes. In agreement with Ref. [7], for weak noise 
the model displays long-range ordered motion (up to the 
actual system size L), that disappears in a continuous 
manner by increasing rj. 

As L — > 00, the numerical results indicate the presence 
of a kinetic phase transition described by 



4>(rf) ~ < V Vc(g,L) 




\/9 

J for 7] < n c {g,L) ^ ^ 
for n > r) c (g, L) 



where i] c (g,L) is the critical noise amplitude that sepa- 
rates the ordered and disordered phases. For a given set 
of 4>{rj), the exponent (3 and t} c {q, L) is determined by se- 
lecting the values providing the best fit to the Ansatz (5), 
i.e., yielding the maximal scaling regime. The numerical 
results are consistent with (5), since the scaling regime is 
increased for larger N (Fig. 2b) and the estimated values 
°f Vc(q, L) converge to a non-zero n c (g, 00) value as 



Vc(Q,L) -r] c {e,oo) ~ N~ 



< 



(6) 



where £ is approximately equal to 0.25 (Fig. 2c). This 
calculation yields (for g = 2) j3 = 0.42 ± 0.03, which 
is definitely different from the the mean-field value 1/2, 
and consistent with the value reported in [7] obtained for 
smaller systems using definition (i) for S(i). 

Next we discuss the role of density. In Fig. 3a, (j>(rf) is 
plotted for L = 100 and various values of g. One can ob- 
serve that the long-range ordered phase is present for any 
g, but for a fixed value of 77, cj> vanishes with decreasing g. 
These 4>{rf) functions parameterized by various g collapse 
to a "universal" function (j){x) by rescaling n with n c (g), 



4>{v,q) = 4>{v/vc{q)), 



(7) 



where <p(x) <~ (1 — x) 13 for x < 1, and <j)(x) « for x > 1, 
and rjc(g) is determined as the value which minimizes 



A(*) = 



J dx (j>(xz, g) — <fi(x)j 



(8) 



Here the A < 1 cutoff is chosen to exclude the noisy 
and rounded (due to finite-size effects) region around rj c . 
In our calculations we used the value A = 0.9, and we 
determined <fr(x) in a self-consistent manner by averaging 
over the already rescaled <j>(rj/rj c (g)) functions. We show 
the result of the data collapse in Fig. 3b. 

This procedure (together with the finite-size analysis 
for a given g) also yields the position of the "critical line" 
rj c (g) in the n — g parameter space. According to our 
numerical results, 



Vc(g) 



(9) 
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holds with k = 0.45 ± 0.05 (see Fig. 3c). Apart from 
the numerical uncertainties, k = 1/2, in agreement with 
recent theoretical results [24]. 

The critical line (9) is qualitatively different from that 
of the diluted ferromagnets, since here the critical density 
at i] — > (corresponding to the percolation threshold for 
diluted ferromagnets, see, e.g., [22]) is vanishing, 



wl(At) = «0 2 ) At - (4>)l t ), 



(13) 



where the internal and external brackets denote averages 
calculated over a time interval [t,t + At] and the en- 
tire stationary data set, respectively. A close relation 
between w^,(At) and the power-spectrum S^,(lo) of the 
order parameter (j>(t) can be derived from the self-afhnc 
properties of the signal [25]: w<f,(At) ~ (At) a is equiva- 



lim Q c {ri) 



0. 



To see this, let us imagine a system of size L consist- 
ing of two particles only. Due to the finite size of S(i), 
for almost all initial conditions the trajectories of these 
particles will get close enough to each other to establish 
interaction. As the noise is negligible in (1), the cluster 
formed by the particles will not break apart, resulting in 
(j) « 1 for any finite g as rj — > 0. 

The behavior of the model in the g — > oo limit is still 
not clear. By definition rj < 2ir holds, so (9) obviously 
cannot describe the system in this limit. Thus t] c (q) ei- 
ther approaches 2ir or a non-trivial r) c (co) < 2ir value. 

Finally we note that Eq.(7) also implies that the ex- 
ponent /?', defined as (f> ~ (g — g c )@ for g > g c (see Ref. 
[7]), must be equal to /?, since 



(10) lent to S<p(ui) ~ uj , and for the exponents A = 1 + 2a 
holds (for 3 > A > 1). 

The numerically obtained results (Fig. 5a) show that 
at the critical line the fluctuations of the order pa- 
rameter are characterized by the correlation exponent 
q = 0.6 ± 0.05 up to a characteristic correlation time 
r. This behavior means that in the steady state the 
'condensation' and 'evaporation' processes (when parti- 
cles join or leave the dominant cluster, respectively) are 
correlated [26]. For t > t the the system shows even 
stronger correlations, as in this regime w^(At) ~ log At 
(see Fig. 5b) . This behavior is probably related to the pe- 
riodic boundary conditions, as r is comparable to L/vo, 
the time needed for a particle to cross the entire system. 



<f>(v,Qc(v) +e) 



n 



l - 



1 dr)c 
r\ dg 



(11) 



where rj c (g) denotes the inverse function of Q c {rf) as 
V = VciQciv)]- Indeed, the results of the simulations 
performed with rj = const and various L and p yield 
/?' = 0.4 ± 0.05, which is consistent with (3 = (3' . In this 
case the larger uncertainty is due to the increased noise 
at low densities. 



III. FLUCTUATIONS 

As a further analogy with equilibrium phase transi- 
tions, we note that the fluctuations of the order param- 
eter also increase on approaching the critical line. To 
study this, we calculate for various control parameters 
the standard deviation of the total momentum, defined 
as a 2 = ((f) 2 ) — {<fi) 2 , where the averages are taken over 
the stationary data set obtained from the simulations. 
In Fig. 4a we plot a vs the rescaled noise amplitude 
x = rj/rj c (g) for various densities and L = 100. The tails 
of the curves are symmetric, and decay as power-laws 
with an exponent 7 close to 2 (see Fig. 4b) 



IV. EXTERNAL FIELD 

' The presence of the long-range correlated phase in two 
dimensions (i.e. the breakdown of the Mcrmin- Wagner 
theorem [23]) is a striking consequence of the non- 
equilibrium nature of the XY model. As the fluctuation- 
response theorem kT\ = Na 2 must hold for Hamiltonian 
systems only, it is interesting to check its applicability 
when the equation of motion cannot be derived from a 
Hamiltonian. 

A natural way to introduce an external field in the 
model is by changing rule (2) for 




(14) 



a{x) 



II 



(12) 



where e is an arbitrary unit vector and the parameter h 
controls the strength of the perturbation. 

In the numerical studies we set the density of the par- 
ticles to q = 2, and study </> as a function of 77 and h. 
Typical results for systems of L = 100 are plotted in 
Fig. 6. The curves parameterized by various 77 intersect 
the h = axis at the same values found in Sec. II, so 
for h = the original model is recovered. In general, 
(j> is increased for increasing h, and the results can be 
summarized by means of critical exponents <5, 7,, 7^ and 
susceptibility x(v) defined in a manner similar to the case 
of classical magnets: 



In Sec. IV we compare this result to direct susceptibility 
measurements, when an external field is also applied. 

We studied the time correlations of the fluctuations by 
calculating the expected value of the rms deviation in a 
time interval At 



xiv) 



lim 



4>( V , h) - 0fa h = o) 

h 



h 1/d for 77 > Tj c , 



(15) 



(16) 
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and 
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-7', 

for 7] < 7] c 



(17) 



We distinguish the critical exponents 7* and 7^ from the 
exponent 7 defined by the singularity of (7(77), since in 
this case 7 ^ 7* ^ 7I (as will be demonstrated later). 

To determine 5, we plot the obtained (f>(r]) curves on 
a double logarithmic plot (see Fig. 7a). We see that for 
large values of <f> the system saturates (</> > <Asat ~ 0.1), 
while for low values of <f> the finite-size noise dominates 
((f) < </>noisc « 0.03, for L = 100). Thus, we expect that 
the relation 



<j>( V ,h)=x(ri)h 1/S 



(18) 



holds for S [(/'noise, 0sat] only. As the boundaries of this 
interval do not depend on rj until rj > r/ c , Eq. (18) can 
be verified by collapsing the rescaled 4>*{i1jh)/x(v) data 
points onto a single power-law (Fig. 7b). Similarly to the 
determination of rj c (g), x(v) an d S can be calculated in a 
self-consistent manner. This procedure yields S = 1.1 a 1 
and x(v) decaying with 7* w 4. Note, that 7* cannot be 
obtained by this method as in the ordered phase <j) > 4>s&t ■ 
To check the stability of these results, we also per- 
formed measurements on larger systems, where </5> n oisc is 
reduced. Indeed, the scaling regime for (18) increases 
and the exponent is consistent with our former estimate 
(see Fig. 7c) 

We can also calculate xiv) a Pplying its definition (17) 
by sampling a series of x{v)i = [<f>{v, M _ HVi K)\l {hi ~ 
fi'i) for various hi, h\ <C 1, and (as S w 1) assuming x{f]) = 
(Xi{v))i- The result is shown on Fig. 8. The 77 > r\ c tail is 
consistent with our former estimate from applying (18), 
while for 77 < r\ c x(j]) decays with 7* « 1. 



V. CONCLUSIONS 

We have demonstrated that this far-from-equilibrium 
system of self-propelled particles can be described using 
the framework of classical critical phenomena, but shows 
surprising new features when compared to the analogous 
equilibrium systems. The velocity vq provides a control 
parameter which switches between an equilibrium XY 
type model (v = 0), and another universality class of 
dynamic, non-equilibrium models characterized by a non- 
vanishing value of f • 

Indeed, for v = we can observe Kosterlitz-Thouless 
vortices in the system, which turned out to be unstable 
for any nonzero vq we investigated: in Fig. 9 we plot 
snapshots of the system. These pictures demonstrate how 
the vortices disappear and give way to long-range order. 
We performed control simulations with various vq in the 
range [0.01,0.3], and the results presented seem to be 
robust against changing the value of vq in this range. 
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FIG. 1. Scematic illustration of the model. The particles 
move off-lattice on a plane and interact with other particles 
located in the local surrounding, which can be either a circle 
or 9 neighboring cells in an underlying lattice. We plot these 
interaction areas for particle A with a solid and dashed line, 
respectively. 

FIG. 2. (a) The average momentum of the system in the 
steady state vs the noise amplitude r\ for g = 2 and four dif- 
ferent system sizes [(o) N = 800, L = 20; (+) N = 3200, 
L = 40; (□) N = 20000, L = 100 and (x) N = 10 5 , L = 223]. 
(b) The order present at small n disappears in a continu- 
ous manner reminiscent of second order phase transitions: 
4> ~ {(r)c(L) - v)/Vc(L)] 13 = {A V /r, c f, with p = 0.42, dif- 
ferent from the mean-field value 1/2 (solid line), (c) The 
estimated rj c (L) values converge to a non-zero ?7 c (oo) limit as 
L -1 / 2 , indicating the presence of the ordered phase even as 
N — > oo. The data represent time averages of long (> 10 5 
MCS) simulations. 

FIG. 3. (a) The average momentum of the system in the 
steady state vs n for L = 100 and three different densities 
[(□) g = 4, (O) g = 2 and (+) g = 0.5]. (b) The <j>(r]) 
functions parameterized by various g can be collapsed onto 
a single curve 4>[q / rj c {g)] = <f>(r),g). (c) The critical line in 
the rj — g phase space is a power-law in the examined regime: 
T)c{g) ~ Q K with n w 0.45 (solid line) for a system of size 
L = 100. 

FIG. 4. (a) The rms deviation of the order parameter a in 
the steady state, for L — 100, and two independent data sets 
(O, +), each averaged over runs with g — 0.3, 0.6, 0.9, 1.2, 1.5 
and 2 . (b) The divergence is symmetric and its tail decays 
as o~(x) ~ |1 — x\~ 2 (solid line), where x denotes the rescaled 
noise amplitude r)/r) c {g). 



FIG. 5. The expected value of the standard deviation of 
the order parameter in a time interval of length At on double 
logarithmic (a) and log- linear (b) plots. For At < r ~ L/vo 
the system shows long-range correlations characterized by 
w r j,(At) ~ At a where a ~ 0.6. For comparison, the dot- 
ted line shows the uncorrelated case (a — 1/2). For At > r 
the correlations are even stronger: w$(At) ~ log(At). The 
data displayed is an average over two independent runs, each 
was performed at n — 2.73, g — 2 and N = 10 5 . The duration 
of the simulation was 70 000 Monte-Carlo steps. 

FIG. 6. <f> vs the amplitude of the applied external field h, 
for g = 2, L = 100 and jj = 0.3,0.9,1.5,2.1,2.7,3.3,3.5 and 
4.5 (from top to bottom). 



FIG. 7. (a) The <j>(h) functions for g = 2, L = 100 and 
tj — 3.9, 4.2, 4.5, 4.8, 5.2 (from top to bottom) on a double-log- 
arithmic plot. <f> ~ h 1 / 6 holds for a limited range of <j> only, 
as for cf> > 0.1 the system saturates, and for (j> < 0.02 the 
finite-size noise dominates, (b) For 0.02 < cf> < 0.1 the curves 
for various n can be collapsed onto a single power-law with 
an exponent 1/ifss 0.9. (c) For larger systems the finite-size 
noise is reduced and the scaling regime is enlarged, in agree- 
ment with the (f> ~ h 1/s Ansatz [(O) g = 2,L = 100, n = 4.2; 
(+) e = 2,L = 223, 7? = 4.2]. 

FIG. 8. The susceptibility of the system obtained with two 
different methods: (□ with error bars) applying definition 
(15), and (O) from data collapse in the </> sa t > 4> > 0noise 
regime, which is a more precise procedure for n > n c . The 
data suggest that in this case the divergence is asymmetric: 
ji ~ 1 (solid line) and 7, « 4 (dotted line). 

FIG. 9. Snapshots of the time development of a system 
with N = 4000, L = 40 and v = 0.01 at 50 (a), 100 (b), 
400 (c) and 3000 (d) Monte-Carlo steps. First the behavior 
is reminiscent of the equilibrium XY model, where the long 
range order is missing since vortices are present in the system. 
However the vortices are unstable, and finally a self-organized 
long-range order develops. 
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